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Protein degradation via ubiquitination is a major proteolytic mechanism in cells. Once 
a protein is destined for degradation, it is tagged by multiple ubiquitin (Ub) molecules. 
The synthesized polyubiquitin chains can be recognized by the 26S proteosome where 
proteins are degraded. These chains form through multiple ubiquitination cycles that 
are similar to multi-site phosphorylation cycles. As kinases and phosphatases, two 
opposing enzymes (E3 ligases and deubiquitinases DUBs) catalyze (de)ubiquitination 
cycles. Although multi-ubiquitination cycles are fundamental mechanisms of controlling 
protein concentrations within a cell, their dynamics have never been explored. Here, we 
fill this knowledge gap. We show that under permissive physiological conditions, the 
formation of polyubiquitin chain of length greater than two and subsequent degradation 
of the ubiquitinated protein, which is balanced by protein synthesis, can display bistable, 
switch-like responses. Interestingly, the occurrence of bistability becomes pronounced, 
as the chain grows, giving rise to "all-or-none" regulation at the protein levels. We 
give predictions of protein distributions under bistable regime awaiting experimental 
verification. Importantly, we show for the first time that sustained oscillations can robustly 
arise in the process of formation of ubiquitin chain, largely due to the degradation 
of the target protein. This new feature is opposite to the properties of multi-site 
phosphorylation cycles, which are incapable of generating oscillation if the total abundance 
of interconverted protein forms is conserved. We derive structural and kinetic constraints 
for the emergence of oscillations, indicating that a competition between different 
substrate forms and the E3 and DUB is critical for oscillation. Our work provides the first 
detailed elucidation of the dynamical features brought about by different molecular setups 
of the polyubiquitin chain assembly process responsible for protein degradation. 
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INTRODUCTION 

Regulated degradation of proteins has pivotal roles in many cel- 
lular processes including cell cycle progression, transcription, 
signal transduction, differentiation, and apoptosis (King et al., 
1996; Hershko and Ciechanover, 1998; Nguyen et al, 201 1, 2013). 
Ubiquitination via the ubiquitin-proteosome system (UPS) is the 
primary degradation pathway through which the stability and 
levels of cellular proteins are controlled with high specificity. Its 
malfunctioning leads to a number of human diseases (Lonser 
et al., 2003; Guerriero and Brodsky, 2012; Sano and Reed, 2013). 
This is exemplified by the quality control endoplasmic reticulum 
(ER) -associated protein degradation (ERAD) pathway that tar- 
gets misfolded or inappropriately accumulated proteins of the ER 
for ubiquitination and subsequent degradation by the 26S pro- 
teosome. Aberrant upregulation of ER stress induced ERAD can 
lead to apoptosis and ERAD dysfunction has been implicated in 
diseases (Guerriero and Brodsky, 2012; Sano and Reed, 2013). In 
yeast, supplying cells solely with the mutated form of ubiquitin 
(Lys48 to Arg) without the wild-type ubiquitin leads to cell-cycle 
arrest and a lethal phenotype (Finley et al, 1994). Disruption of 



the ubiquitination machinery, particularly the ligating enzymes, 
have also been implicated in the pathogenesis of several dis- 
orders including the Von Hippel-Lindau syndrome, Angelman 
syndrome and the autosomal-recessive growth retardation disor- 
der 3-M syndrome (Lonser et al., 2003; Huber et al., 2005; Greer 
et al., 2010). Recent large-scale quantification of mammalian gene 
expression has revealed significant variability in the stability of 
thousands of cellular proteins with half-life ranging from a few 
minutes to a few months (Schwanhausser et al., 2011), further 
highlighting the versatility of the UPS pathway. 

Monoubiquitination, the covalent tagging of a single ubiq- 
uitin molecule to the target substrate protein, is a multi-step 
process that involves a cascade of three essential enzymes: an 
El ubiquitin- activating enzyme, an E2 ubiquitin conjugating 
enzyme, and an E3 ubiquitin ligase (Figure lA). The C-terminus 
of ubiquitin first forms a thioester bond with the catalytic systeine 
of the El in an ATP-dependent manner before being trans- 
ferred from El to the catalytic systeine of the E2. Subsequently, 
the E3 ligase binds both the ubiquitin-charged E2 and the sub- 
strate to catalyze transfer of the ubiquitin to a lysine residue 
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FIGURE 1 I Biological scliematic diagrams of ubiquitination. 

(A) Scliematic diagram illustrating the sequential reaction steps of ubiquitin 
activation and transferring involving a cascade of E1/E2/E3 enzymes, 
leading to polyubiquitination of the target protein. (B) Simplified schematic 
diagram displaying the formation of two major types of polyubiquitin chain 
via repeated ubiquitination cycles: Lys48 and Lys63-linked chains, with 
differential functional consequences. 



on the substrate to form an isopeptide bond, resulting in sub- 
strate monoubiquitination (Hershko and Ciechanover, 1998). 
Importantly, ubiquitin contains seven lysine residues which can 
be utilized by the appropriate E2/E3 pair to catalyze further 
cycles of ubiquitination, assembling various types of polyubiqui- 
tin chains. Polyubiquitination through ubiquitin Lys48 linkages 
generally targets the substrate protein for degradation, and is 
the principal focus of this work, while other linkages such as 
Lys63-linked chains are associated with non-proteolytic functions 
including kinase activation, signal transduction, and endocyto- 
sis, as illustrated in Figure IB (Ye and Rape, 2009). Similarly to 
phosphorylation cycles, the action of the E3 ligase is reversed by 
a class of enzymes called deubiquitinases (DUBs) which catalyze 
the reverse deubiquitinating reaction. 

A number of recent theoretical investigations have been car- 
ried out to explore the dynamic properties of signaling networks 
where ubiquitination and UPS-mediated degradation play major 
roles in their signal regulation (Nguyen et al., 2011, 2013). These 
studies have revealed that ubiquitination can bring about intri- 
cate dynamic behaviors including bistable switches, oscUlatory 
and excitable responses to the system of interest. In these pre- 
vious studies the assembly of polyubiquitin chains on target 
proteins were always considered as a single step. However, the 
chain assembly is in fact a dynamic and multi-reaction process. 
Here, we provide for the first time a detailed dynamic inves- 
tigation of the assembly of polyubiquitin chain on a substrate 
protein mediated by multiple cycles of ubiquitination and sub- 
sequent degradation. Employing kinetic modeling, simulations 
and bifurcation analysis, we demonstrate that the formation 
of polyubiquitin chains is a highly dynamic process capable of 
generating complex temporal and steady-state behaviors. Using 
parametric and bifurcation analysis, we derive regions of dis- 
tinct dynamics of the system. Further, we examine the effect of 
the ubiquitin chain length on the degradation dynamics and the 



steady-state distribution of the targeted protein's concentration. 
Applying the first-passage theory, we provide intuitive insights 
into the connection between protein concentration and average 
protein lifetime. Our results provide novel mechanistic insights 
into the intricate operation of the polyubiquitin chain assembly 
machinery. 

RESULTS 

A CORE MODEL OF POLYUBIQUITIN CHAIN ASSEMBLY 
Development of a core kinetic model 

Polyubiquitin chain assembly is an essential marking mecha- 
nism for the recognition of a protein substrate destined for 
destruction. Previous studies have shown that the assembled 
chain must consist of at least four ubiquitin molecules for effi- 
cient proteasomal targeting (Thrower et al., 2000). Two compet- 
ing models for the chain extension currently exist. The chain 
is either buQt in a stepwise manner with sequential addition 
of one ubiquitin monomer after another, or following an "en- 
bloc" transfer wherein the ubiquitin chain is pre-assembled on 
a E2 enzyme before being transferred to the substrate in a sin- 
gle step (Li et al., 2007; Ravid and Hochstrasser, 2007; Ye and 
Rape, 2009). Furthermore, a combined model where chains are 
formed following both mechanisms in a mixed manner can also 
exist. 

For simplicity, we fist consider a core model of polyubiqui- 
tin chain assembly consisting of only two ubiquitination cycles 
while longer cycles are examined in later sections. A schematic 
scheme describing the mechanistic interactions of this core model 
is given in Figure 2A. Here E3 ligase is assumed as a main catalyz- 
ing enzyme, regulating the dynamics of the chain formation in 
this model, while neglecting the presence of El and E2 enzymes. 
This assumption is supported by data suggesting that E2s, such as 
Cdc34, are mostly pre-charged with ubiquitin under steady-state 
condition (lin et al, 2007). Importantly, the final ubiquitinated 
substrate with the longest assembled chain (S-ub2) is assumed to 
be subject to rapid degradation by the 26S proteosome, whereas 
the unmodified substrate or intermediate ubiquitinated substrate 
forms are under a basal, substantially lower degradation rate 
{kct >> kh). The substrate degradation is balanced by a constant 
synthesis of the unmodified substrate S. We further assume that 
ubiquitin is in abundance (and not a limiting factor), there- 
fore its concentration is considered a constant parameter in the 
core model. Based on this mechanistic reaction scheme, we con- 
structed the model using ordinary differential equations (ODEs). 
We describe the (de)ubiquitination catalyzed by the E3 ligase and 
its opposing DUB enzyme at elementary step levels that foUow 
mass-action kinetics which. Model reactions and ODEs are given 
in detail in Tables 1A,B. 

Since polyubiquitin chain assembly can follow combined 
extension modes via stepwise or en-bloc transfer (Raman and 
Harper, 2009), chains with four ubiquitins can be formed from 
two ubiquitination cycles in several ways; e.g., by stepwise trans- 
fer followed by an en-bloc transfer of a three-ubiquitin chain, 
or two sequential en-bloc transfers of two-ubiquitin chains. This 
implies that our two-cycle core model described above is relevant 
for analysing the dynamics of chain assembly with length of four 
or even more ubiquitin monomers. 
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FIGURE 2 I Bistable, "all-or-none" responses in the core model. 

(A) Mass-action kinetic sclieme of a two-cycle chain formation system 
regulated by E3/DUB. The reactions and ordinary differential equations 
(ODEs) of the core model are given in Tables 1A,B. (B) Quasi steady state 
(QSS) curves intersection on the [S] vs. [S-ub2l plane determines the 
system steady states: 3 intersection points for bistability vs. 1 intersection 
point for monostability. (C,D) Bistable response and hysteresis curves of 
the system: dependence of steady-state levels of [S] and [S-ub2] on 
increasing abundance of the E3 ligase. The blue and red solid lines indicate 
stable and unstable steady states, respectively. The vertical dashed lines 
illustrate the systems steady states corresponding to three cases of E3 
total: monostability (Mi and M2) and bistability (Bi_3). Parameters values 
and units used for plotting are given in Table 1A, 4th column. 



Key experimental observations of ubiquitination Idnetics 

Experiments aimed at obtaining kinetic insights and/or differenti- 
ating between competing models of ubiquitination are challeng- 
ing due to the fast speed of ubiquitin-transfer reactions (Pierce 
et al., 2009; Raman and Harper, 2009). Nevertheless, a num- 
ber of key experimental observations have been made which 
narrow our choice of the kinetic parameter values for the core 
model. Using the cyclin-dependent kinase inhibitor Sicl as a 
substrate and its cognate E2/E3 pair (the cuUin-RING ubiqui- 
tin ligase SCF'^'''^* and the ubiquitin-conjugating enzyme Cdc34), 
Petroski and Deshaies have examined the Lys48 polyubiquitina- 
tion process (Petroski and Deshaies, 2005). They revealed that 
Sicl polyubiquitination occurs in two distinct stages: addition of 
the first ubiquitin which is slow and rate limiting, followed by the 
attachment of subsequent ubiquitins at a much faster rate. As a 
result, Sicl containing a single ubiquitin molecule (Sicl-ub) was 
ubiquitinated up to 5-fold faster than unmodified Sicl (Petroski 
and Deshaies, 2005). This indicates that assuming the same DUB 
kinetics for the two ubiquitination cycles in the core model, the 
flux of the first E3-catalyzed ubiquitination cycle would be sig- 
nificantly lower than that of the second cycle. Moreover since 
SCF was found to increase the Vmax of the di-ubiquitin elon- 
gation reaction (Petroski and Deshaies, 2005), we can assume 
that the differential ubiquitination rate was primarily due to the 
lower catalytic rate (fccat) of the first cycle. More recent data 
lend further support to this observation (Pierce et al., 2009). 
Employing the same pair of E2/E3 but using cyclin E (CycE) 



and P-Catenin (P-Cat) as substrates. Pierce et al. performed mbc- 
ing reactions on a millisecond timescale in a quench-flow device 
which allowed them to subsequently visualize the different sub- 
strate forms (Pierce et al., 2009). Fitting their data to a theoretical 
chain building model based on sequential transfers of ubiquitins, 
the authors could estimate the individual transfer and dissoci- 
ation rates for each intermediate in the generation of polyu- 
biquitinated CycE/p-Cat. Consistent with the previous results, 
they found that for both substrates once they are monoubiqui- 
tinated, the rate of chain elongation dramatically increases before 
slowing down again as the chain grows further. Interestingly, 
the substrate dissociation rate values are similar regardless of 
ubiquitin-chain lengths (Pierce et al., 2009). Overall, these 
data have provided useful qualitative knowledge for our model 
parameterization. 

To further quantitatively parameterise the model, we estimated 
the core model kinetic parameters by equating the lifetime of the 
substrate between our model and Pierce et al's. Derived values 
for the E3-related parameters are given in Table 2. In the absence 
of well documented kinetic information on DUBs, we assume 
that DUBs operate on a similar timescale as the E3 ligases and 
DUB-related kinetic parameters have the same values for different 
ubiquitination cycles. 

BISTABLE "ALL-OR-NONE " SWITCHES IN PROTEIN DEGRADATION 
Emergence of bistable switches in the presence of degradation 

We have previously reported the existence of bistable switches 
in a dual phosphorylation-dephosphorylation cycle with a dis- 
tributive mechanism for the catalyzing kinase and phosphatase 
(Markevich et al., 2004). Bistability is a phenomenon in which 
a dynamic system switches between two distinct stable steady 
states but cannot rest in an intermediate state. Bistable behav- 
iors have been observed in many biological systems and shown 
to be involved in cell-fate decision making and cell-cycle con- 
trol (Bagowski and FerreU, 2001; Xiong and FerreU, 2003; Wang 
et al, 2009). A fundamental difference that distinguishes multi- 
ple phosphorylations from polyubiquitination (such as on Lys48) 
is the absence of rapid ubiquitination-triggered degradation. 
We thus explore if protein degradation would still give rise to 
bistable switches in the synthesis of polyubiquitin chains and 
if so, how would bistability be influenced by the length of the 
chain. 

Analysis of the core model using parameter values derived 
in the previous sections show that abrupt, bistable switches 
can also emerge in the presence of rapid degradation of S-uba 
(Figure 2). The occurrence of bistability can be illustrated by plot- 
ting the nuUclines (i.e., the solution curves for the ODE right 
hand sides wherethe time derivative is set to zero) on the sys- 
tems phase plane. The intersection of the nuUclines gives the 
equilibrium points of the ODE system. Strictly speaking, nuU- 
clines describe two-dimensional systems. In high-dimensional 
systems including our model, bistability can similarly be demon- 
strated by plotting the quasi steady-state (QSS) curves for two 
chosen variables (e.g., two substrate forms, S and S-uba) on one 
plane where the system's steady state also correspond to inter- 
section points of the curves. Figure 2B shows that depending 
on the total level of the E3 ligase, there can either be one or 
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Table 1A | Reactions and reaction rates of the core, double-cycle model in Figure 2A. 



No Reactions Reaction rates Parameter values for Figure 2 Parameter values for Figure 4 



1 


S + E3 S-E3 


Vl = 


<rfi ■ S ■ E3 - kn 


- S-E3 


^fi 


= 0.06, kn 


= 0.1 


kn 


= 0.002, kn - 


= 0.7 


2 


S-E3 E3+S-ub 


V2 = 


<rci ■ S-E3 




fei 


= 0.2 




fccl 


= 0.2 




3 


S-ub + DUB Sub-DUB 


V3 = 


<r,2- S-ub- DUB 


- k,2- Sub-DUB 


kn 


= 0.01 , kr2 


= 0.1 


kf2 


= 0.01 , k,2 = 


0.006 


4 


Sub-DUB DUB+S 


V4 = 


<rc2Sub-DUB 




kc2 


= 0.2 




kc2 


= 0.25 




5 


S-ub + E3 « Sub-E3 


V5 = 


ki3- S-ub- E3 - 


kr3- Sub-E3 


kn 


= 0.01, kn 


= 0.1 


kn 


= 0.01, kn = 


0.1 


6 


Sub-E3^E3+S-ub2 


V6 = 


<rc3-Sub-E3 




kc3 


= 4 




kc3 


= 4 




7 


S-ub2 + DUB Sub2-DUB 


v? = 


<rf4-S-ub2-DUB 


- /Cr4-Sub2-DUB 


ki4 


= 0.01 , kr4 


= 0.1 


kn 


= 0.06, fcr4 = 


0.043 


8 


Sub2-DUB ^ DUB+S-ub 


V8 = 


<rc4-Sub2-DUB 




kc4 


= 0.2 






= 0.054 




9 


0 ^ S 


Vg = 


ks-S 




ks-- 


= 0.08 




ks-- 


= 0.08 




10 


S-ub2 ^ 0 


VlO = 


- fcd-S-ub2 




kd- 


= 0.01 




kd-. 


= 0.01 




11 


0 


Vll = 


-.kfS 




kb- 


= 0.0001 




kt- 


= 0.0001 




12 


S-ub^ 0 


Vl2 = 


= /ff,-S-ub 

















The first- (dissociation kr, catalytic kc, degradation k^, kt) and second-order (association kf) rate constants are expressed in s ^ and nM ' s '. Synthesis rate is 
expressed in nMs'' ' 



Table 1B | Ordinary differential equations of the double-cycle model. 



Left-hand sides Right-hand sides Initial concentrations (nM) 



d[Sl/dt 


Vg - 


Vl -1- 


V4 - 


Vll 


10 


d[S-E31/dt 


Vi - 


V2 






0 


d[Sub-E31/dt 


V5 - 


V6 






0 


d[S-ub]/dt 


V2 - 


vs - 


Vg -1- 


Vg - Vl2 


0 


d[S-ub2l/clt 


V6 - 


V7 - 


VlO 




0 


d[Sub2-DUB]/dt 


v? - 


va 






0 


d[Sub-DUB]/dt 


V3 - 


V4 






0 


d[E31/dt 


-Vl 


-h V2 


- V5 


-I-V6 


100 


d[DUB]/dt 


-V3 


-h V4 


- V7 


-l-Vg 


100 



The reaction rates are given in Table lA. 



three intersection points between the QSS curves which corre- 
spond to monostable or bistable responses, respectively. At an 
intermediate level of E3 ligase, the QSS curves (solid lines) inter- 
sect at three points Bi_3 resulting in bistability. The lower or 
higher total E3 shifts the QSS curve for S down (dashed) or up 
(dashdotted) on the plane, while leaving the S-ub2 QSS curve 
intact, resulting in just one intersection point (Mi or M2) and 
monostability. 

As a consequence of bistability, hysteresis is observed for most 
system species, illustrated in Figures 2C,D for the unmodified 
substrate S and final ubiquitinated form S-ub2 in response to 
increasing E3 level (Figures 2C,D). The specific E3 total values 
selected for plotting in Figure 2B and the respective intersec- 
tion points are shown on the hysteresis curves. As the E3 ligase 
abundance gradually increases from Mi, the system starts at a 
high steady-state S (low steady-state S-ub2) and follows the high 
branch into the bistable regime before abruptly switching to a 
low S value (high S-ub2) at an E3 threshold value character- 
ized by a saddle-node bifurcation. On other other hand, if the 
total E3 begins at a high level (M2) and gradually decreases, 
the system will traverse the low steady-state S (high S-ub2) 
branch into the bistable regime before abruptly switching to 



the high branch at a second, lower E3 threshold value defined 
by another saddle-node bifurcation point. Moreover, with E3 
total within the bistable regime, either the low or high stable 
branch (Bi or B3) can be reached depending on the system's 
history. 

Substrate production and degradation rates affect bistabiiity 

To examine how the model parameters directly involved in pro- 
tein level regulation affect bistability, we analyse the system 
behavior in response to change in the rate of production (fcj) 
and degradation (fcj) of the protein substrate. Summing up all 
the substrate forms (ubiquitinated and unubiquitinated) avail- 
able in the model, we also observed hysteresis for the total 
S against fcj and fcrf. However, opposite changes occur when 
either kg or kd increase. (Figures 3A,B). Bistability was found 
over a wide range of the production rate, with the upper 
bound exceeding the normal physiological range, suggesting 
that sufficiently high production rate is necessary for bistable 
behavior (Figure 3A). On the other hand, although bistabil- 
ity persists over a large range of degradation rate, exceedingly 
strong degradation abolishes bistability (Figure 3B). The QSS 
plots for the corresponding various dynamics scenarios that 
are illustrated in Figures 3A,B are also shown in Figures 3C,D, 
where changes in the parameters fcj and fcj shift both QSS 
curves. 

ft should be noted that in the absence of any degradation 
(even with a very small basal rate for the unmodified substrate 
S, kh = 0), under the specific conditions of low E3 ligase activity, 
the concentration of S can increase unlimitedly without reaching 
a steady state within a physiologically relevant timeframe. This 
is because the rate of S production would outweigh that of its 
maximal depletion rate. Although such case cannot be ruled out 
theoretically, it is highly unlikely to occur in vivo. Moreover, we 
found that if only the final substrate form, S-ub2 is degraded, 
a double (de)ubiquitination cycle cannot display bistability, as 
the final ubiquitinated form S-ub2 becomes independent of other 
variables at the steady state and always exhibits a single steady- 
state value. Thus, the inclusion of a basal degradation rate (which 
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Table 2 | Model reactions, equations and kinetic parameter values for simplified models with varying chain length from 2 to 5 in Figures 5, 8. 



Reactions 

E3 catalyzed ubiquitination S(-ubn) binding to E3 ligase 

dissociation of S(-ubn)-E3 complex 
Catalysis 



Substrate S synthesis rate, kg 

degradation rate, kf, 
degradation rate, kfj 



Unit 1st cycle 2nd cycle 3rd cycle 4th cycle 5th cycle 



nM-i s-i 0.01 0.01 0.01 0.01 0.01 

s"'' 0.1 0.1 0.1 0.1 0.1 

S"'' 0.2 4 1 0.6 0.3 



0.01 0.01 0.01 0.01 0.01 

0.1 0.1 0.1 0.1 0.1 

0.2 0.2 0.2 0.2 0.2 



nM s-i 0.05 _ _ _ _ 

s-i 0.0002 0.0002 0.0002 0.0002 

s-i - - - - 0.02 



DUB catalyzed de-ubiquitination S(-ubn) binding to DUB deubiquitinase nM~ s~ 

dissociation of S(-ubn)-DUB complex s~^ 
Catalysis s~^ 



The ODEs for the models of length 1-5 can be constructed straighforwardly from the schematic diagram In Figure 5A in a similar way as described in Tables 1A,B 

for model of length 2. 
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FIGURE 3 I Dependence of bistability on model parameters. (A,B) 

Bistable response and hysteresis curves showing dependence of the 
steady-state levels of [S] and |S-ub2] on the rate of production (ks) 
and degradation (k^) of the substrate protein. The blue and red 
solid lines indicate stable and unstable steady states, respectively. 
The vertical dashed lines illustrate the systems steady states 
corresponding to three values of (A) ks and (B) k^ displayed in 
(CD). (C,D) Corresponding QSS curves intersection determines the 
system steady states: three intersection points for bistability vs. 1 
intersection point for monostability. The remaining parameters values 
and units used for plotting are given in Table 1A, 4th column. 



may be mediated by processes other than the UPS, e.g., the 
lysosome or other cellular proteases) for the intermediate forms 
of the substrate is necessary for the system to achieve bistable 
behavior. 

Bistability for non-proteoiytic polyubiquitination 

Non-proteolytic polyubiquitination, exemplified by the Lys63- 
linked ubiquitin chain, has been known not to target the 
substrate for rapid degradation but instead direct other cellu- 
lar functions. In this case, the chain assembly process would 



be an un-leaking, so-called "closed" system which resembles 
multi-site phosphorylation. Similarly to Lys48 polyubiquiti- 
nation, bistability can also emerge in this system [data not 
shown, as this system is mathematically equivalent to a multisite 
phosphorylation system previously shown to exhibit bistability 
(Markevich et al., 2004)]. 

EMERGENCE OF SUSTAINED OSCILLATIONS DUE TO RAPID 
DEGRADATION 

Emergence of sustained oscillations in the core model 

Despite the absence of any explicit negative feedback in our core 
model of chain assembly, we hypothesized that the existence of 
bistable behavior coupled with protein degradation which serves 
as a sink can bring about sustained oscillatory dynamics. We thus 
set out to examine this hypothesis. Comprehensive exploration 
of the parameter space, using systematic parameter sampling 
and visualization of the temporal dynamics, has enabled us to 
identify parameter domains where the system displays oscilla- 
tions. Further analysis showed that robust oscillations can emerge 
within physiologically plausible parameter ranges (Figure 4). 

Figure 4E shows the temporal oscillatory dynamics of S and 
S-ub2 at different levels of the total E3 ligase. The E3 level strongly 
influences the shape of oscillations, producing wave-like (with 
period in minutes) to distinct pulse-like pattern (period in hours) 
as in the case of S-ub2 concentration as the level of E3 decreases. 
Quantifying the oscillation amplitude and period revealed that 
both increase exponentially with decreasing E3 concentration 
until a threshold level is reached where the oscillations terminate 
abruptly (Figures 4C,D). To further characterize the oscillation 
dynamics, we carried out both local and global bifurcation anal- 
yses using Mathematica and XPPaut (Wolfram Research, 2010). 
A one-parameter bifurcation diagram in Figure 4A shows a sta- 
ble limit cycle arising from the Hopf bifurcation accompanying 
an unstable fixed-point (dashed). This limit cycle suddenly dis- 
appears as the E3 concentration decreases below a specific value, 
denoted by the left red dot in Figure 4A. Furthermore, a two- 
parameter bifurcation plot can partition the E3-DUB plane into 
distinct dynamic regions, revealing a large oscillation regime 
indicative of the robust emergence of this behavior (Figure 4B). 
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FIGURE 4 I Sustained oscillations in the core model. (A) Bifurcation 
diagram of steady-state lS-ub2l concentration (blue line) vs. total E3 level 
showing existence of oscillatory dynamics in the core model. Dashed blue 
indicates the region of oscillations, and the red dashed lines indicate the 
minimum and maximum value of the oscillatory dynamics. (B) Two-parameter 
partition of system dynamics projected onto the E3-DUB plane. (C) 
Dependence of the oscillation amplitude (in log scale) on the total E3 level, 
blue lines indicate the minimum and maximum value of the oscillatory 



dynamics (as in A). (D) Dependence of the oscillation period (in log scale) on 
the total E3 level. The left and right dashed lines indicate the E3 values for 
which the oscillatory dynamics emerge and terminate. (E) Distinct temporal 
patterns of system species ranging from wave-like to pulse-like oscillatory 
dynamics resulted from different levels of the total E3 ligase. Parameters 
values used for plotting are given in Table 1A, 5th column. The x-marks 
indicate the maximum and minimum amplitude of the oscillations as used to 
generate (A,C). 



It is important to point out that a sufficient degradation rate 
of the polyubiquitinated substrate is an essential condition for the 
observed oscillations. Our analysis shows that lack of degrada- 
tion or even very small degradation rates prevent the oscillations. 
This suggests that oscillatory behavior is an intrinsic property of 
the "open" polyubiquitination-degradation system (i.e., the total 
abundance of interconverted protein forms is not conserved), 
that does not occur in "closed" systems such as multi-site 
phosphorylation which do not feature rapid degradation. 

EFFECTS OF CHAIN LENGTH ON SYSTEMS DYNAMICS 

Although our two-cycle core model has provided novel insights 
into the dynamic behavior of ubiquitination chain extension, 
polyubiquitin chains are often made up of multiple ubiquitin 
monomers with length up to 10 (Pierce et al., 2009). Here, we 
investigate how both bistability and oscillations found in the 
two-cycle model are affected as chain elongation takes place. 

Bistability is enhanced as polyubiquitin chain elongates 

The increase in the number of ubiquitination cycles results in 
an increased non-linearity in the system and consequently in a 
steeper dose response against E3 ligase, which in turn enhances 
switch-like behavior of the system (Figure 5). This effect of the 
chain length on the dose response behavior is analogs to the result 
for a canonical MAPK cascade with multi-site phosphorylation 
cycles: each cascade layer increases steepness of the downstream 
response to external stimulation (Brown et al., 1997; Ferrell, 1997; 



Kholodenko et al., 1997). We find that for a biologically-derived 
parameter set (Table 2) the total substrate concentration indeed 
responds in an increasingly more switch-like fashion to changing 
E3 ligase levels in the presence of a longer ubiquitination chain. 
Moreover, owing to faster ubiquitination rate in the second cycle 
(fccat = 4s^') compared to latter cycles (fccat < Is^'), bistabil- 
ity emerges for ubiquitination chains longer than 2 (Figure 5B). 
The width of the bistable region increases only slightly for 
longer chains, which might indicate that the switch-like reg- 
ulation of the substrate rather than the hysteretic response 
with wide bistable region is functionally relevant for cellular 
physiology. 

Closely associated with steep switch-like response and bistabil- 
ity is the notion of a bimodal distribution of the protein level in a 
clonal population. In our system, each cell differs to some degree 
with respect to E3 and DUB concentration due to divergent gene 
expression profiles. Consequently, the overall population-level 
response consists of a superposition of slightly differing indi- 
vidual dose-response curves, which may give rise to a bimodal 
distribution of the substrate concentration. Such a response on 
the level of a cellular population may occur when the E3 enzyme 
level differs between cells within the ultrasensitive regime of the 
sigmoidal dose response (Birtwistle et al., 2012). Alternatively, 
bimodality may arise when E3 levels vary within the bistable 
region of the hysteretic dose-response (Becskei et al, 2001; Paliwal 
et al., 2007). We illustrate the latter in Figure 5C for the 5-cycle 
ubiquitination chain. Such population-level measurements with 
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FIGURE 5 I Dependence of the dose response on chain length. (A) 

Mass-action kinetic scheme of models consisting of 1 up to 5 ubiquitination 
cycles. Note that the reactions and ODEs for the models of length >3 are 
not explicitly given but can be easily constructed from the kinetic scheme 
following the format outlined in Tables 1A,B for length 2. (B) Dose- 
response curves for the total S concentration show the emergence of 
bistability for chains longer than 2 ubiquitination cycles. Dotted lines 
indicate the bistable region. (C) Steady-state distribution of total substrate 
concentration in all forms for the ubiquitination chain of length 5. Total E3 
and DUB concentrations exhibit cell-to-cell variability modeled by gamma 
distribution with standard deviation equal to 10 nM. The histogram was 
obtained by solving numerically the system of corresponding ODEs for 
pairs of E3/DUB concentrations sampled 10,000 times from this 
distribution. Parameter values used for plotting are given in Table 2. 



single-cell resolution can be performed with flow-cytometry or 
high content screening. 

Relation between protein levels and average protein lifetime 

Stochasticity inherent to all biochemical processes at the level 
of a single cell introduces variability in the protein lifetime. 
The lifetime x of molecules undergoing simple one-step degra- 
dation follows an exponential distribution with the coefficient 
of variation (CV = SD / mean) equal 1 (Li and Qian, 2002). 
A sequential degradation process such as ubiquitination steps 
discussed here may result in a peaked, gamma-like lifetime dis- 
tribution, which becomes narrower with an increasing number 
of steps (Dobrzynski, 2008; Pedraza and Paulsson, 2008; Shin 
et al, 2009). As a consequence the CV decreases too. In the sim- 
plest case of an irreversible reaction sequence with equal rate 
constants for each of the steps, the CV depends only on the num- 
ber of steps N and follows CV = 1/(VN) (Dobrzynski, 2008). 
This narrowing down of the lifetime distribution has been argued 
to bring about precision in biological processes: a much tighter 
regulation may take place when molecules exist in the cell for 
a duration narrowly distributed around the mean. For exam- 
ple, Li and Qian have demonstrated that GTP hydrolysis cycle 
modeled as four reversible steps significantly reduces the vari- 
ance in the lifetime of the GTP-bound state (Li and Qian, 2002). 
Such an improved accuracy is arguably advantageous for pre- 
cise signal amplification to downstream effectors. Due to very 
similar lifetimes (from a narrow distribution), each GTP-bound 
molecule activates an almost equal number of effectors. Signal 



relay is therefore robust every time the receptor signal triggers 
production of few GTP-bound molecules. 

Peaked and narrow lifetime distribution also affects the decay 
pattern of the molecule (Dobrzynski, 2008). The analysis by 
Deneke et al. (2013) of experimentally measured mRNA decay 
time courses revealed deviation from the exponential decay and 
the existence of fast and slow decay patterns. These in turn have 
been attributed to sequential mRNA decapping that gives rise to 
non-exponential lifetime distributions. The effect of such distri- 
butions on the steady-state protein or mRNA variability is subtler, 
however. A seminal result from queuing theory states that for sys- 
tems with exponentially distributed synthesis events, the lifetime 
distribution has no effect and the resulting variability of molecule 
level is Poissonian; such systems are known to be insensitive to 
the degradation ("service") time distribution (Adan and Jacques, 
2001). However, the result may rarely hold in a biological setting 
for the exponential timing of protein or mRNA synthesis often 
occurs in bursts, which strongly affect the stationary distribution 
(Shahrezaei and Swain, 2008). 

We calculate protein lifetime distributions for various ubiq- 
uitination chain lengths shown in Figure 5A using first-passage 
theory (Redner, 2001). We conclude that the distributions are 
peaked but their CVs are close to 1, as in exponential distribu- 
tion, and decrease weakly (less than 1%) with the addition of 
a single ubiquitination cycle (Figure 8). The multi-step degra- 
dation process is therefore close to memoryless and the protein 
lifetime distribution is well approximated by an exponential. As 
a consequence, the theoretical prediction follows that the decay 
pattern of a protein targeted for multi-step ubiquitination will 
undergo the classic exponential decay, for which the protein half- 
life equals Ln(2) x < t > (angle brackets denote the average 
over all molecules of a given protein species within a cell, i.e., the 
average protein lifetime). Pulse-chase technique (Schwanhausser 
et al, 2011) or light activated fluorescent tag such as Dendra2 
(Zhang et al., 2007) can be used to measure decay dynamics. 

Regardless of the time distribution of synthesis and degrada- 
tion events, the average protein lifetime, < t >, can be calculated 
by applying another general result of queuing theory, the Little's 
Law (Elgart et al., 2010). According to the theorem, the aver- 
age lifetime is related to the average total protein level < S > 
and the average protein synthesis rate Xs through the relation 
<S>=Xs<T>.A useful ramification of this relation is that 
any of the three quantities can be calculated by measuring the 
remaining two. More importantly, the rule holds regardless of 
the statistics governing protein production and degradation. In 
addition, the Little's rule implies that the bistability evidenced by 
the hysteretic dose-response in total S concentration shown in 
Figure 5B exists also in the average protein lifetime — the obser- 
vation corroborated by the calculation of lifetime distributions 
shown in Figure 8B, where two lifetime distributions emerge that 
correspond to two stable protein level states. 

BISTABILITY AND OSCILLATIONS PERSIST IN EXTENDED MODELS 

Development of the detailed E1/E2/E3 model 

As discussed earlier, ubiquitination is a multi-step process involv- 
ing not only the ligases but also the El activating enzyme, the E2 
conjugating enzymes and the availability of ubiquitin molecule 
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FIGURE 6 I Oscillations in the extended E1/E2/E3 model. (A) Kinetic 
diagram sliowing tlie detailed model of polyubiquitin chain formation of 
length 2, with explicit incorporation of ubiquitin (Ub) and the El, E2 
enzymes. The reactions and ODEs of the this model are given in 
Tables 3A,B. (B) Dependence of oscillatory patterns on the abundance of 
available ubiquitin ranging from low to high level in the system where Ub is 
conserved. Parameters values and units used for plotting are given in 
Table 3A. 



itself (Figure 1). It is thus important to see whether complex 
behaviors such as bistable and oscillatory responses still persist 
when potential effects exerted by the E1/E2 enzymes and ubiqui- 
tin status are explicitly accounted for. To this end, we extended 
the core model to incorporate the action of El, E2, and ubiqui- 
tin within a detailed mechanistic model of two-cycle ubiquitin 
chain assembly, whose kinetic scheme is shown in Figure 6A. The 
model scheme follows largely the biological sequence of reactions 
described at a mechanistic level in Figure 1. A number of key 
assumptions were made. 

First, we assumed that ubiquitin (Ub) becomes activated when 
bound to El, and then it is loaded onto the E2 enzyme. The 
E2~Ub complex then binds to the pre-assembled S-E3 com- 
plex and transfers Ub to the substrate, releasing free E2. The 
E3 ligase can also dissociate from the complex yielding the free 
substrate in ubiquitinated form. The whole cycle repeats when 
the second Ub is conjugated to the mono-ubiquitinated sub- 
strate or when further Ub molecules are attached. Similar to the 
core model, the substrate is constantly produced and subjects 
to a low basal degradation rate while the end-product S-ub2 is 
rapidly degraded by the 26S proteosome. Furthermore, ubiqui- 
tinated substrates are continuously deubiquitinated by the DUB 
resulting in Ub recycling. Implied by the name, ubiquitin is 
ubiquitously expressed in most eukaryotic tissues and is a sta- 
ble protein due to its unique globular structure (Shabek and 
Ciechanover, 2010). The in vivo pool of ubiquitin is therefore 
likely to be maintained at homeostasis under normal conditions 
by the dynamic equilibrium between multiple processes of Ub 
degradation and Ub synthesis. However, such pool can be altered 
in specific pathological contexts (Shabek and Ciechanover, 2010). 
In light of literature evidence, we thus consider two relevant sce- 
narios: (1) when Ub is conserved in the system and (2) when 



Ub can be degraded together with the ubiquitinated substrate 
(Shabek et al., 2009) in which case the Ub pool is balanced 
by a constant synthesis of Ub (Figure 6A). The model reac- 
tions, ODEs and parameter values for both scenarios are given 
in Tables 3A,B. 

Oscillations and bistability in the extended model 

We found that sustained oscillations persist in this more realistic 
model of chain building. Interestingly, the abundance of ubiqui- 
tin was found to strongly control the existence of oscillations and 
the oscillatory pattern (Figure 6B). In the case where Ub level 
is conserved, oscillations only arise when sufficient Ub is avail- 
able. However, there appears to be no upper threshold of the Ub 
pool for the existence of oscillations, suggesting more abundant 
ubiquitin facilitates this dynamics. Moreover, higher Ub level is 
associated with lower oscillation period (Figure 6B). Similar to 
what was observed in the core model system, we also found that 
oscillations can only occur within defined (bounded) ranges of 
the E3 and DUB total levels. Similar observations were made for 
the El and E2 abundances, implying that appropriate levels of 
these key enzymes are necessary for the oscillation to arise. 

In the second scenario where Ub can be degraded by the 26S 
proteosome together with the substrate, the total Ub availabil- 
ity is a dynamic variable. Nevertheless, model analysis showed 
that to maintain oscillations, the Ub level should also be prop- 
erly supplied as too low Ub abundance eliminates oscillation. 
Figure 7 A shows an initial temporal oscillation of model species 
which subsequently vanishes when the level of the depleting 
Ub pool becomes lower than the threshold required for oscilla- 
tions (Figure 7B). Increasing the Ub production rate can rescue 
oscillations by compensating for the level of lost Ub. 

Furthermore, bistable behavior also exhibits in the two- 
cycle extended models in both cases when ubiquitin is con- 
served or being continually turned over. Consistent with the 
results from the core model, the bistable regime occurs over 
a defined range of E3 and DUB total and is also sensitive 
to the substrate production and degradation rates. We further 
extended the model to include up to four ubiquitination cycles. 
Importantly, having more ubiquitination cycles as a result of 
chain elongation does not appear to compromise the existence 
of both bistability and oscillations. Taken together, the results 
here indicate that oscillations and bistability are still inher- 
ent to the polyubiquitin chain formation process even when 
the effects of other enzymes than the E3 ligase are explicitly 
incorporated. 

DISCUSSION AND CONCLUSION 

COMPETITION OF DIFFERENT SUBSTRATE FORMS FOR THE 

CATALYZING ENZYMES IS CRITICAL FOR COMPLEX DYNAMIC 

BEHAVIORS 

In the models of ubiquitin chain assembly presented above, sub- 
strate reactions catalyzed by the regulating enzymes E3 ligase 
and DUB are modeled explicitly using a mass-action descrip- 
tion at the elementary step level. Under this kinetic formula- 
tion, the different forms of the target substrate compete for 
the catalyzing enzymes. For instance, the increased E3 portion 
bound in the first ubiquitination cycle wiU result in less E3 
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Table 3A | Reactions and reaction rates of the two-cycle extended model in Figure 6A. 



No. 


Reactions 


Reaction rates 


Parameter values 


1 


S + E3 S-E3 


vi = /tfvS.ES - <rn-S-E3 


kfi = 0.01, kn = 0.1 


2 


S-E3 + E2-ub^ E2 + 


V2 = kav S-E3.E2-ub 


fcai = 0.00003 




Sub-E3 






3 


Sub-E3 ^ E3 + S-ub 


V3 = kci ■Sub-E3 


fcci = 0.2 


4 


S-ub + DUB Sub-DUB 


V4 = kf2- S-ub.DUB - /fr2-Sub-DUB 


kf2 = 1.5, kr2 = 0.15 


5 


Sub-DUB^ DUB + S + 


V5 = /rc2-Sub-DUB 


fcc2 = 1 


6 


ub 

S-ub + E3 ^ Sub-E3 


Vg = /Cf3- S-ub. E3 


% = 1 


7 


Sub-E3 + E2-ub^ E2 + 


V7 = /ra3-Sub-E3.E2-ub 


fca3 = 0.01 




Sub2-E3 






8 


Sub2-E3 ^ E3 + S-ub2 


vg = /rc3-Sub2-E3 


fcc3=4 


9 


S-ub2 + DUB 


Vg = /tf4.S-ub2-DUB - /(rr4-Sub2-DUB 


tf4 = 1.5, k,4 = 0.15 




Sub2-DUB 






10 


Sub2-DUB DUB + 


Vio = /Cc4-Sub2-DUB 


kcA = 0.04 




S-ub + ub 






11 


El -h ub El-ub 


Vii = fcfs.EI .ub - /fr5-E1-ub 


tf5 = 0.1, <rr5 = 0.06 


12 


El-ub-h E2^ El + E2-ub 


Vi2 = /r6-E1-ub-E2 


kfi = 0.004 


13 


0^8 


Vi3 = ks-S 


ks = 0.06 


14 


0 ^ ub (*) 


Vi4 = ks2-S 


ks2 = 0.001 /ts2 = 0 (*) 


15 


S-ub2 -> ub -h ub 


Vi5 = kci-S-ub2 


fcd = 0.05 




S-ub2 ^ 0 (*) 






16 


S^ 0 


Vie = kjj-S 


kt = 0.0001 


17 


S-ub^ 0 


Vi7 = fct,.S-ub 





77ie first-order (dissociation kr, catalytic kc, degradation k^, kt) and second-order (association kfj rate constants are expressed in ' and nl\/t ' ' . Syntliesis rate 
is expressed innlVI s^K 

"Note that for tfie modei wliere Ub is degraded witti substrate, reaction 14 is absent (ks2 = 0), and in reaction 15 Ub is not recycled. 



Table 3B | Ordinary differential equations of the two-cycle extended 
model. 



Left-hand sides Right-hand sides 



Initial concentrations (nIVI) 



dlSl/dt 


Vl3 - 


-Vl6 ■ 


■1-V5 - 


- Vl 


10 


cllS-E31/dt 


Vi - 


V2 






0 


d[Sub-E31/dt 


V2 - 


V3 






0 


d[S-ubl/dt 


V3 - 


V4 - 


V6 -h 


Vio - Vl7 


0 


d[Sub2-E3]/dt 


V7 - 


V8 






0 


d[S-ub2l/dt 


V8 - 


Vg - 


Vl5 




0 


d[Sub2-DUBl/dt 


Vg - 


Vio 






0 


d[Sub-DUBl/dt 


V4 - 


V5 






0 


d[E3]/dt 


-Vl 


-h V3 - 


-Ve- 


fVs 


100 


dIDUBl/dt 


-V4 


+ V5- 


-Vg- 


1-Vio 


100 


d[E1]/dt 


-Vll 


+ Vi; 


? 




1000 


d[E1-ubl/dt 


Vll - 


- Vi2 






0 


d[E2]/dt 


-Vl2 


-h V2 


+ V7 




1000 


d[E2-ubl/dt 


-Vl2 


- V2 


- V7 




0 


d[ub]/dt (*) 


Vl4 - 


-Vll - 


fVg- 


)-Vio -|-2Vi5 


10,000 



The reaction rates are given in Table 3A. 

"Note that for the model where ubiquitin is conserved, equation (») is not 
included. 

available for the second cycle. To examine if such competi- 
tion should be a prerequisite for the emergence of oscillations 
and/or bistability, we modeled instead the (de)attachment of 
a single Ub in the core model with simple non-competitive 
Michaelis-Menten (MM) formulation, where the reaction rate 



implicitly contains the concentration of the catalysing enzyme 
[Vmax = fccat* (Enzyme)]. Simple non-competitive MM descrip- 
tion of both E3 and DUB-regulated reactions failed to lead 
to oscillations or bistable switches. However, using a form 
of competitive MM description based on total quasi-state 
approximation described in our previous work (Markevich 
et al., 2004), bistability and oscillations can both be found. 
These results together indicate that enzyme competition is 
critical for the existence of oscillations and bistability, as 
such competitions coupled with proper level of non-linearity 
would endow the system with effect of positive or nega- 
tive feedback which underlie bistable response and oscillations, 
respectively. 

SUPPORTING EXPERIMENTAL EVIDENCE OF SWITCH-LIKE 
UBIQUITINATION 

The theoretical predictions of potentially bistable and oscilla- 
tory behavior await future experimental validation. The field 
of ubiquitination, ubiquitin-related signaling and ubiquitin- 
associated detection techniques is maturing with rapidly acceler- 
ating pace. The studies of the epidermal growth factor receptor 
(EGFR)-activated signal transduction (where ubiquitination is 
strongly involved in signaling regulation (Nguyen et al, 2013), 
have shown that in response to graded levels of the EGF 
ligand, EGFR ubiquitination follows a switch-like, threshold- 
characterized response (Sigismund et al., 2013). This observations 
were verified by the measurements of ubiquitinated EGFR using 
both ELISA and stable isotope labeling with amino acids in 
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cell culture (SILAC)-based mass-spectrometry. Importantly, the 
switch-like dose-response of the ubiquitinated EGFR correlates 
exactly with the EGFR internalization mediated by non-clathrin 
endocytosis, a primary mechanism responsible for EGFR degra- 
dation (Sigismund et al., 2013). 

In vitro experimental design based on the model analysis 
results could also be the next step in confirming the predictions 
about the dynamics of the systems components. A reconsti- 
tuted in vitro system with purified forms of relevant proteins 
(i.e., E3 ligase and DUBs) provides flexibility where one can 
explore wide ranges of enzyme concentrations. If compartmen- 
tal localization is required, some of these involved proteins may 
be embedded into a phospholipid membrane bUayer or lipo- 
somes, which also can be utilized to facilitate complexes forma- 
tion and enhance reactions rates. Potential temporal oscillatory 
behavior can be explored by adding the substrate protein fol- 
lowed by the addition of ubiquitin, the E1/E2 enzymes, E3 lig- 
ase and DUB to the reaction medium. At periodically selected 
time points, aliquots are taken, and the ubiquitinated level of 
the substrate can be measured by immunoblotting using rele- 
vant antibodies. Specific antibodies raised against common types 
of linkage including Lys48 and Lys63 are now available. More 
direct in vivo approaches such as imaging techniques using 
microscopy-based binding assay can also be exploited for high 
temporal resolution measurement [105], which have been previ- 
ously used to detect oscillations in GTPases network (Carlin et al., 
2011). Rapid development of quantitative mass-spectrometry 
based methods for measuring ubiquitinated protein level in bulk 
would be another viable option (Sigismund et al., 2013). In 
these in vivo setups, one may choose to inhibit the 26S pro- 
teosome to prevent UPS-mediated degradation, and therefore 
better identify the ubiquitination signals. On the other hand, 
detection of switches can be done by similar measurement tech- 
niques in response to increasing titration of a dose component, 
for example the regulating E3 ligase. Single-cell based approaches 
such as flow-cytometry or high content screening can serve an 
ideal platform for detecting the bimodal distribution of pro- 
tein levels, indicative of the existence of switches or bistable 
responses. 

IMPLICATIONS OF SWITCH-LIKE PROTEIN DEGRADATION REGULATION 

There are multiple accounts of switch-like regulation that 
have been demonstrated experimentally for biological systems. 



Cascaded signaling networks exhibit steep dose responses that 
can translate analog input signal into all-or-none response of the 
transcription factor (Brown et al., 1997). Likewise, cooperative 
binding of transcription factors can give rise to steep relation 
between the factor concentration and the promoter activity (Rossi 
et al., 2000; Veitia, 2003). Here, we suggest an additional layer 
of digital regulation at the level of protein degradation due to 
a sequence of ubiquitination cycles. Our theoretical considera- 
tions demonstrate that for proteins degraded by way of multi- 
step polyubiquitination, steep or even hysteretic dose responses 
may arise when changing the parameters, such as the synthe- 
sis and degradation rates or the total ligase and deubiquitinase 
concentrations. Hence, a number of properties beneficial for 
cellular physiology can be argued. For example, even without 
pronounced steepness in the dose response at the signaling level 
or without strong cooperativity at the stage of transcription 
initiation, the protein level can stiU be regulated in a switch- 
like manner in response to intra- or extra-cellular cues. Such a 
switch-like dependence may also facilitate a filter to temporal 
fluctuations in the rate of protein synthesis introduced by tran- 
scriptional and translational bursting or pausing. Furthermore, 
resistance to fluctuations extends to changes in the concentra- 
tion of the ligase and deubiquitinase (Figure 5B). Due to the 
steep sigmoidal relation between total substrate S and E3/DUB, 
large variations in the latter translate to digital responses in 
S. Interestingly, for the biologically-relevant parameters we find 
that the switch-like response is enhanced in the presence of 
longer ubiquitination chains. Notably, only proteins tagged with 
at least four ubiquitins are targeted by the proteasome, which 
strengthens our conjecture that switch-like regulation of protein 
degradation resulting from bistability may indeed be function- 
ally utilized by cells to facilitate cellular decisions. Since there are 
less than 1000 different types of E3 ligase and even far less for 
E2 and El enzymes (Ye and Rape, 2009), a common E3 ligase 
is likely to be shared by multiple target substrates. In this case, 
a degradation regulation system enabled by multiple switch-like 
dose-response curves with different E3 threshold values for each 
substrate could be a possible mechanism to generate substrate 
specificity. 

CONCLUDING REMARK 

Since the discovery of protein ubiquitination more than three 
decades ago, subsequent work has revolutionized our per- 
ception of its role in signaling networks. As a principal 
mechanism for specific cellular protein degradation control- 
ling many key cellular processes, polyubiquitination exhibits 
overwhelming diversity and complexity whose functional prop- 
erties only begin to be unraveled. Our work contributes 
to this endeavor by providing a comprehensive investigation 
into the dynamic features of the process of ubiquitin chain 
formation. 

MATERIALS AND METHODS 
CONSTRUCTION OF THE KINETIC MODELS 

To construct the kinetic models analyzed in section Results 
which includes: (1) the simplified, double-cycle ubiquitination 
model used in section A core model of polyubiquitin chain 
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FIGURE 8 I Lifetime distributions of the substrate protein S. (A-C) 

Dashed lines show the numerically calculated distributions of S for the 
models depicted in Figure 5. Solid lines show the exponential 
distributions with the same average lifetime as the corresponding lifetime 
distribution. Note that the exponential function is a very good 



approximation and discrepancies occur only for short lifetimes. Bistability 
of the system with chains longer than 2 results in two lifetime 
distributions that correspond to two stable states of the system, shown 
in (B). For clarity, we show the curve for four ubiquitination cycles 
only in (C). 



assembly, Bistable "all-or-none" switches in protein degradation, 
and Emergence of sustained oscillations due to rapid degra- 
dation (2) the simplified models with increasing chain length 
(up to 5 cycles) used in section Effects of chain length on sys- 
tems dynamics and (3) and the detailed E1/E2/E3 model used 
in section Bistability and oscillations persist in extended mod- 
els, we employed a deterministic modeling approach based on 
ordinary differential equations. The reaction rates were for- 
mulated at the elementary step level using mass-action kinetic 
law. Detailed reactions, reaction rates and the ODEs as well 
as parameter values for the models are given in the fol- 
lowing tables. Specifically, the simplified, double-cycle model 
can be reconstructed from Tables 1A,B; the simplified models 
with varying chain length can be reconstructed from Table 2 
and the detailed E1/E2/E3 model can be constructed from 
Tables 3A,B. 

ESTIMATION OF KINETIC PARAMETERS FOR SEQUENTIAL 
UBIQUITINATION CYCLES 

The lack of experimentally measured kinetic data remains a 
general challenge for computational modeling. For our mod- 
els of sequential ubiquitination cycles, the choice of kinetic 
parameters has been guided by key experimental observa- 
tions wherever available or typical values for protein associa- 
tion/dissociation and enzymatic reaction rates, as explained in 
section Key experimental observations of ubiquitination kinet- 
ics. In particular, we estimated the kinetic parameters for the 
simplified models of various lengths, given in Table 2, based 
on in vitro measurement by Pierce et al. (2009). We assume 
that the rate of catalysis, fcc(n) in our model, corresponds 
to the ubiquitin transfer rate by Cdc34-SCF E3 ligase. We 
set complex dissociation rate, kr(„) in our model, to 0.1 s^', 
which is of the order of the dissociation rate /Coff measured by 
Pierce et al. 

SIMULATION AND BIFURCATION ANALYSIS PROCEDURE 

Models implementation, temporal simulations, and calculation 
of the QSS curves were conducted in Mathematica. QSS curves 
were calculated as previously described (Nguyen et al., 2011). 



Bifurcation diagrams (including saddle-node and Hopf bifur- 
cation) and the dependencies of steady states on parameters 
were calculated by numerical continuation algorithms using the 
software XPPAUT (www.math.pitt.edu/~bard/xpp/xpp.html). 

CALCULATION OF PROTEIN LIFETIME DISTRIBUTIONS USING 
FIRST-PASSAGE THEORY 

After synthesis, a single substrate molecule S undergoes a 
sequence of biochemical reactions such as formation of an 
intermediate complex with E3 ligase or deubiquitnase, and addi- 
tion of ubiquitin tags (Figure 5A). In our model, the lifetime of 
S is the time between the synthesis and the degradation through 
either of the channels with rate constant ki, or /cj. These reactions 
are of stochastic nature due to molecular fluctuations within the 
cell, hence the lifetime of S will follow a distribution rather than 
a single fixed number. In order to calculate that distribution we 
resort to first-passage time theory (FPT) (Redner, 2001). 

First, we note that the ODEs that describe the dynamics of our 
system (e.g.. Tables IB or 3B) correspond to the set of chemical 
master equations for subsequent states of substrate S modifica- 
tions. The rates of transitions between the states are described 
by reaction rate constants multiplied by steady-state concentra- 
tions of relevant chemical species. For instance, the transition 
from state S to enzyme-bound S-E3 takes place with the rate kfi x 
[S—E3]. By taking the Laplace transform of both sides of the origi- 
nal time-dependent system, we obtain a set of algebraic equations, 
which we solve for all the state variables: S, Sub, S-E3, etc. The life- 
time of S is the first time the single molecule S reaches any of the 
degradation states (signified by dots in Figure 5A) from the ini- 
tial state [S] = 1. Using the example of a two-cycle ubiquitination, 
the FPT is obtained as a flux from all the states that may directly 
lead to the degradation. In the Laplace-transformed s-space, 
flux is: 

/ (s) = h ([S] (s) + [Sub] (5)) + kdSub2](s) 

By Taylor-expanding /(sj we can calculate moments of the FPT 
(or protein lifetime) and obtain the mean or variance. The inverse 
Laplace transform will yield the time-dependent function which 
has the interpretation of the probability density function. We 
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plot it for various lengths of the ubiquitination sequence in 
Figure 8. 
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